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We study the synchronization phenomena in a system of globally coupled oscillators with time 
delay in the coupling. The self-consistency equations for the order parameter are derived, which 
depend explicitly on the amount of delay. Analysis of these equations reveals that the system in 
general exhibits discontinuous transitions in addition to the usual continuous transition, between the 
incoherent state and a multitude of coherent states with different synchronization frequencies. In 
particular, the phase diagram is obtained on the plane of the coupling strength and the delay time, 
and ubiquity of multistability as well as suppression of the synchronization frequency is manifested. 
Numerical simulations are also performed to give consistent results. 
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d . I. INTRODUCTION 

a 

■ When a large population of limit cycle oscillators with slightly different natural frequencies are coupled, they often 
come to oscillate with an identical frequency. Such collective synchronization phenomena have been observed in 
Q various oscillatory systems in physics, biology, chemistry, and other sciences Q-Qj, attracting much interest in recent 
years p[-pj]|. In the Kuramoto model for those oscillator systems, oscillators are coupled with each other via the 
interaction which depends on the phase difference between each pair || . It describes the emergence of phase coherence 
' with the increase of the coupling strength, elucidating interesting connection between the collective synchronization 
. and a phase transition. 

Here, as in the usual dynamics of many-particle systems, sufficient attention has not been paid to the effects of time 
retardation in the oscillator system. In biological systems such as pacemaker cells and neurons, however, temporal 
delay is natural and the finite time interval required for the information transmission between two elements may 
be important [^|^]. Time delay in the interaction may modify drastically dynamic behavior of the system, such as 
stability and ergodicity [^2[ . In some types of a system of coupled oscillators, retarded interactions have been found 
| to result in multistability and suppression of the collective frequency |l3|-|l6[]. In a system of two coupled oscillators, 
it has been found that the time delay induces a multitude of synchronized solutions. Namely, in the system with 
finite time delay, more than one stable solution are possible at given coupling strength. Among those, the most 
stable solution is the one with the largest synchronization frequency, as shown via the linear stability analysis [ jl3| . 
Similar behaviors including frequency suppression and multistability have been observed in the neural network model, 
where peripheral oscillators with identical natural frequencies are coupled only with a central oscillator by forward 
and backward connections with time delay ||l4| . The two-dimensional system of identical oscillators with time-delayed 
nearest- neighbor coupling has also been considered to reveal similar frequency suppression [jl5| . The system of non- 
identical oscillators with delayed interactions has been studied recently |17|| : In the case of a Lorenzian distribution 
of natural frequencies with a nonzero mean, the stability boundary of the incoherent state has been obtained and 
coexistence of one or more coherent and/or incoherent states has been observed in appropriate regions [p"7| . However, 
detailed behaviors such as frequency suppression and emergence of different coherent states have not been addressed 
fully. 

This paper investigates in detail the effects of time delay in the interaction on collective synchronization of coupled 
oscillators with different natural frequencies. For this purpose, we derive the self-consistency equations for the order 
parameter and examine how the characteristic features of the collective synchronization change due to time delay. This 
reveals a multitude of coherent states with nonzero synchronization frequencies, each separated from the incoherent 
state by a discontinuous transition. In particular, we show that the system with a nonzero average frequency can be 
reduced to the system with the vanishing average natural frequency, which allows us to focus on the latter system. 
As in the system without delay, there exists the critical coupling strength, at which the system undergoes the usual 
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continuous transition from the incoherent state to the coherent state displaying collective synchronization (with zero 
synchronization frequency). In addition, at higher values of the coupling strength (beyond the critical value), coherent 
states with larger synchronization frequencies also appear via discontinuous transitions. Thus coherent states with 
different synchronization frequencies in general coexist in the appropriate regions, leading to multistability. The 
synchronization frequency of the oscillators in a coherent state is observed to decrease with the delay time, which is 
similar to the result of other systems with time delay fl4|-[l6f| . 

There are five sections in this paper: Section II presents the system of globally coupled oscillators with time delay, 
as a generalization of the Kuramoto model. The stationary probability distribution for the system is obtained, and 
the self-consistency equations for the order parameter are derived. Section III is devoted to the analysis of the self- 
consistency equations, which reveals the characteristic behavior of the system as the coupling strength or the delay 
time is varied. In particular, the phase diagram is obtained on the plane of the coupling strength and the delay time, 
and ubiquity of multistability as well as suppression of the synchronization frequency is demonstrated. Numerical 
simulations arc also performed and the results, which are in general consistent with the analytical ones, are presented 
in Sec. IV. Finally, Sec. V summarizes the main results, while some details of the calculations are presented in 
Appendices A and B. 



II. SYSTEM OF COUPLED OSCILLATORS WITH TIME DELAY 



The set of equations of motion for N coupled oscillators, each described by its phase 0, (i = 1, 2, • • • , N), is given 

by 

= w< - MMt) ~ Mt-r)}, (1) 

i 

where the prime restricts the summation such that j ' ^ i. The first term on the right-hand side represents the natural 
frequency of the ith oscillator, which is distributed according to the distribution function g{uj). Here giyi) is assumed 
to be smooth and symmetric about a->o, which may be taken to be zero without loss of generality (see below), and 
also to be concave at u) = 0, i.e., g"(0) < 0. The second term denotes the global coupling of strength K/N between 
oscillators, with time delay, indicating that each oscillator interacts with other oscillators only after the retardation 
time t. Without time delay, Eq. (|l|) exactly reduces to the Kuramoto model. 

In order to describe collective synchronization of such an N oscillator system, we define the complex order parameter, 
whose amplitude represents the degree of synchronization, to be 



JV 

i 



I^ e %=Ae w . (2) 



N 

3=1 

Here it is convenient to introduce new variables tpi defined by ipi = 4>i — ilt, where ft is a constant. Note the existence 
of physical invariance due to the rotational symmetry of the total system. In terms of the new variables, Eq. ([!]) reads 

& = -^2J sil #*(*) ~ ^j(* _t ) + fir L ( 3 ) 

3 

where uj i = ten — il. Multiplying Eq. (|J) by e~ jnt , we also obtain the corresponding order parameter for the new 
variables 



1 N 

*^E e * = Ae ^ ( 4 ) 

where 6 = 6 — fit. Incidentally, the order parameter defined in Eq. (|J) allows us to reduce Eq. (||) into a single 
decoupled equation with time delay 

ipi = Qi - KAsm(ipi - 9o), (5) 

where #o = 6 — fir. Although A, f2, and #o depend on the delay time r, they are assumed to be independent of time 
t, which is possible due to the symmetry Considering the relation between the old order parameter and the new one 
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* = *e 4i ", (6) 

we understand that the collective synchronization can be described in terms of a giant oscillator rotating with the 
frequency Q which is in general nonzero. For finite delay time, there exists a multitude of synchronized solutions with 
nonzero values of O; this is in contrast to system without delay, where the rotational symmetry of the system allows 
us to set f2 = 0. 

Instead of Eq. (|J), which may be regarded as a Langevin equation without noise, one may resort to the corresponding 
Fokker-Planck equation for the probability distribution P(ip,t) at zero temperature jTq |: 

-P(?M) = ^[KAsin^ -9 a )-Co}P(^,t). (7) 



The order parameter given by Eq. (|4|) then obtains the form 

N 

X 



N ^ 

i=i 

so 

dCo g{uj+Cl) (e'^t-rjo, (8) 

DO 

where Co — co — O has been noted, and the average of e 1 ^ is to be taken over the distribution P{ip,t—T) of Eq. (]?]) 
with given Co: (e**)t_ r . = fjf dtp P(ip,t-T) e** . 

In the stationary state, we take the average over the stationary distribution P' ' (-0; Co) of Eq. (^). With the 
stationary solution 

( 6[tp — - sirr^ /if A)] for |w|<ifA 

p (0) OA;^)^ ^Co 2 -{ka) 2 . (9) 

^ — n: r~~r — — t- — — t otherwise, 

I 27r|£-A:Asin(V>-0o)| 

it is easy to compute the average: 

(e**)a>= d^P^^-Co) & *+ 
Jo 

( i(Co/KA)-i v / (Co/KA) 2 -1, w > if A 

= e l9 ° ^ i(Co/KA) + y T^ (g/jf A) 2 , -if A < a) < if A (10) 

i i(Co/KA) + iy/(Co/KA) 2 - 1, cD < —if A. 

It is thus natural to divide the system into two groups: one satisfying \Co\ < if A, which is called the synchronization 
group, and the other \Co\ > if A, the desynchronization group. Accordingly, we write 

A = A s + Ad, (11) 

where A s /d is the contribution from the synchronization/desynchronization group to the order parameter. The 
contribution from the synchronization group is given by 



A s e mr = KAj dx g{n+KAx) [y/l - x 2 + ix], 



(12) 

where x = Co/KA. Separating Eq. ( |12[ ) into the real and imaginary parts, we obtain the two coupled nonlinear 
equations 

A s cos fir = if A J dx g(Q,+KAx) \Jl-x 2 , 

A s sin ftr = if A y dxg(n+K Ax) x. (13) 
Similarly, the desynchronization group leads to the equation 
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dx g(Q+KAx) (x + \/x 2 - ij 



4KA 



(14) 



Ad cos 57r = 0, 
A d sin fir = KA 



dxg(tt+KAx) (x + \J x 2 - l) 



(15) 



Note that, unlike the Kuramoto model, the imaginary parts of Eqs. (|T^) and (|l4|) do not vanish, which arises from 
the fact that the nonzero collective frequency due to time delay breaks the symmetry of the integration interval of the 
distribution g{uj). It is obvious in Eq. ( |l5| ) that the contribution from the desynchronization group vanishes in the 
absence of time delay (r = 0) . Recalling that in the presence of time delay the total order parameter is given by the 
sum of the two contributions, one from the synchronizatio n g roup and the other from the desynchronization group, 
we finally obtain the self-consistency equations from Eqs. ( |l3|) and (|T^): 

A cos Qt = KA J dx g(Q+KAx) \/l-x 2 , 

J dxg(n+KAx) x + J dxg(n+KAx) (x - y/a 

/-l 
dxg(Q+KAx) (x + \J x 2 - l" 
-oo ^ ^ 



A sin fir = KA 



(16) 



When the average natural frequency is not zero (loq ^ 0), we define the variables tpi = <f>i — (Q + to^jt, and obtain 
exactly Eq. (|l^) except for f2 replaced by fi = Q + loq. For example, the first equation is given by 



pi 

A cos Qf = KA J dx g(Q.+KAx) \/l - x 2 , 



(17) 



where f is the delay time and the distribution g(uj) is symmetric about uj = loq. Since the distribution g(uj) = g(uj+ujQ) 
is symmetric about lo — 0, we rewrite the above equation in the form 



Acos[(n+uj Q )f} = KA J dxg{n+KAx) VT 



(18) 



which, with the identification (fi+ajo)T = fir, just reproduces the first equation in Eq. ([Tq). Accordingly, the behavior 
of the system with wo ^ 0, which has been considered in Ref. JT7| , can be obtained from that of the system with 
luq = via appropriate rescaling of parameters. 



III. ANALYSIS OF THE SELF-CONSISTENCY EQUATIONS 

The non- vanishing imaginary part of the self-consistency equation given by Eq. (|l6|), arising from time delay, leads 
to a variety of behaviors which are not displayed by the system without delay. In this section, we solve the two coupled 
equations in Eq. ( |l6| ) to obtain the synchronization frequency Q and the order parameter A. We take the Gaussian 

distribution with zero mean and unit variance for the natural frequencies: <?(w) = (27r) _1 / 2 e _w / 2 , and first compute 
numerically the synchronization frequency and the order parameter for various values of the coupling strength K and 
the delay time r. 

Figure 1 exhibits the dependence of the synchronization frequency on K and r, which manifests multistability. 
At small values of r, only the nontrivial solution (A ^ 0) with fi = appears for K > K c (w 1.596), as in the 
system without delay. For large r, on the other hand, solutions with fi ^ also emerge as K is increased further. 
In Fig. 1(a), where the synchronization frequency is plotted as a function of the delay time r at K = 10, it is also 
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observed that the synchronization frequency is suppressed as time delay is increased. This is expected since the delay 
tends to disturb synchronization fl4|| . In Fig. 1(b), we plot the synchronization frequency as a function of the coupling 
strength at r = 5. It shows that at given values of r the synchronization frequency Cl depends rather weakly on K 
after synchronization sets in. Among those solutions at given coupling strength K, the most stable solution is the 
one with the largest value of Cl |L3| although the basin of attraction in general shrinks with Cl. 

The phase boundaries separating the coherent states (A ^ 0) with various synchronization frequencies from the 
incoherent state (A = 0) are shown in Fig. 2, where data have been taken with the step width St = 0.06. Below the 
lowest boundary, which is the straight line K — K c 1.596, only the incoherent state exists; above it, the coherent 
state with O = also exist. Similarly, above each boundary in Fig. 2(a), a new additional coherent state with a 
larger synchronization frequency emerges. Note here that the region of the existence of coherent states constitutes 
two-dimensional (semi-infinite) surfaces in the three-dimensional (K,t,CI) space. Whereas Figs. 1(a) and (b) may be 
regarded as the cross-sections of these surfaces at given values of K and of t, respectively, Fig. 2(a) represents the 
projection of the boundaries of these surfaces onto the K-t plane. In Fig. 2(b) the curves of Fig. 2(a) are redrawn, 
with the horizontal axis rescaled: f in (b) corresponds to CIt/(CI + 3) in (a). In this new scale of the horizontal 
axis, unlike in Fig. 2(a), the boundaries intersect with each other, and only the envelope consisting of the curve 
segments with lowest values of K at given f is displayed in Fig. 2(b). According to the discussion in Sec. II, Fig. 2(b) 
describes the phase boundary below which only the incoherent state exists in the system with loq — 3. Above the 
boundary, the coherent state with the appropriate (nonzero) synchronization frequency, depending on f, appears 
and can coexist with the incoherent state. Note that the lowest boundary (K = K c ss 1.596) in Fig. 2(a) has no 
counterpart in Fig. 2(b) since the zero synchronization frequency (Cl — 0) corresponds to f = 0. Similar boundaries 
have been obtained through numerical simulations for the Lorenzian as well as the delta- function distributions [ ji7| . 
It is of interest to compare Fig. 2(b) with Fig. 4 in Ref. jl7jj, which indicates that the Gaussian distribution leads to 
smoother stability boundaries than the Lorenzian distribution. 

In Fig. 3, the obtained order parameter A is depicted as a function of K at r = 5. Each line describes the transition 
for each synchronization frequency, and the critical coupling strength K c is shown to increase for the transition with 
larger Cl. For example, the leftmost curve, which corresponds to O = 0, gives K c w 1.596, whereas the next one, 
corresponding to O k 1.09, gives K c ss 1.97. Note that as K approaches K c (k, 1.596), the order parameter with 
Cl = decreases continuously to zero, indicating that the leftmost curve describes a continuous transition at K c . 
On the other hand, the rest of the curves with CI > apparently display jumps in the order parameter, indicating 
discontinuous transitions. Accordingly, whereas the lowest boundary in Fig. 2(a) describes a continuous transition, 
the others as well as the boundary in Fig. 2(b) correspond to discontinuous transitions. 

To understand the nature of these transitions analytically, we assume KA <C 1 near the transition to the coherent 
state and expand Eq. (|l^) to the order of (KA) 3 , together with Cl also expanded accordingly: 

Cl^Cla + ^KA + Cl^KA) 2 . (19) 

We investigate two regimes, Cl <C KA and Cl 3> KA, which exhibit phase transitions of different types with each 
other, still taking the Gaussian distribution with zero mean and unit variance for g(uj). 
For Cl <C KA (<§; 1), the self-consistency equation for the order parameter takes the form 

A = ai KA + b^KA) 2 + Cl {KAf + O(KA) 4 , (20) 

where the coefficients a\, b%, and c\ depend on CIq, Cl\, and CI2 defined in Eq. (|l9|). Their specific forms as well as the 
details of the calculation are given in Appendix A. Since the condition Cl <C KA implies CIq <§C 1, we need to consider 

Eq. ( A. 3 ) only for the range —n/2 < tan -1 ^-\/2/7ri7oe n o/ 2 ^ < n/2. It is then obvious that the desired solution of 



Eq. (A. 3) is simply CIq = Cl\ = CI2 = 0, regardless of r. Inserting these values into Eq. (A.£), we obtain the values of 
the coefficients: a\ — 0.626, 61 = 0, and c x = —0.078. Figure 4 illustrates the graphical solution of Eq. (p0|), displaying 
/(A) = [a\K — 1)A + bi (KA) 2 + c\(KA) 3 versus A for b\ — 0. Note that the critical coupling strength is given by 
K c (= ai _1 ) sa 1.595, which indeed agrees perfectly with the numerical value given by the leftmost curve in Fig. 3. 
It is thus concluded that the system displays a continuous phase transition with O = 0, which is consistent with the 
result of the Kuramoto model. 

In the opposite case of O ^ KA (<g; 1), we can still obtain the self-consistency equation for A up to the order of 
(KA) 3 , in a manner similar to that for the previous small-0 case: 

A = a 2 KA + b 2 (KA) 2 + c 2 (KA) 3 + O(KA) 4 , (21) 

where the coefficients a 2 , b 2 , and c 2 again depend on CIq, Cli, and Cl 2 (see Appendix B). In this case, we need to obtain 
larger solutions, considering the regions (n + l/2)7r < tan -1 ^-y/2/7r0 e°o/ 2 ^ < (n + 3/2)7r with non-negative integer 
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n. Interestingly, this in general yields nonzero values of f2o and accordingly, nonzero values of 62, with which Eq. ( |2l| ) 
displays a jump in A at K = — 4c2(&2 2 — ^a^) -1 , thus indicating a discontinuous transition 0]. Such discontinuous 
transitions are ubiquitous for regions with higher values of n. Namely, the system with delay is in general characterized 
by nonzero values of the synchronization frequency together with discontinuous transitions, which is consistent with 
the numerical results displayed in Fig. 3. There the jumps in A displayed by the curves with Q > 0, associated with 
discontinuous transitions, may invalidate the assumption KA <C 1, and the expansion in Eq. is not expected to 
yield quantitatively accurate results. Nevertheless the appearance of such discontinuous transitions has been revealed 
by the above expansion, which is concluded to give a qualitatively correct description of the nature of transitions. 

For more accurate results, we investigate these transition phenomena by examining in detail the behaviors of the 
solutions of the self-consistency equations. Note that A = is always a solution of Eq. for all values of O. To 
seek for other solutions, we divide Eq. (|l^) by A and obtain 



K- 1 cosnr = J dxg(tt+KAx) y/l- x 2 , 

/OO pOO 
dxg(n+KAx)x~ dx g(Q,+KAx) \/ x 2 - 1 
-00 Ji 

L 

dxg(n+KAx) \/x 2 - 1, (22) 

> 

which may be computed numerically. The resulting values of Q versus A are plotted in Figs. 5-7 for r = 5. The solid 
and the dotted lines represent solutions of the first (real part) and the second (imaginary part) equations of Eq. (0), 
respectively. In each figure, the point where the two lines meet with each other provides the synchronization frequency 
and the order parameter A. Figure 5(a) shows the absence of a meeting point for K = 1.58, which implies that 
synchronization does not set in yet. In contrast, the meeting of the solid and dotted lines is obvious for K = 1.60 
shown in (c); (b) reveals a continuous transition (for fl = 0) at K — K c w 1.596, which coincides with the previous 
result. The value of A grows continuously as K is increased beyond K c (see Fig. 6). When K reaches the value 
1.97, as displayed in Fig. 6(b), there emerges via a tangent bifurcation an additional meeting point at finite values of 
A (« 0.08) and Vt {k, 1.09), giving rise to a discontinuous transition, in agreement with the result shown in Fig. 3. As 
the coupling strength is increased further, there appear two meeting points, giving two values of the order parameter 
for the pair of the lines (i.e., with almost the same value of O), as shown in Fig. 6(c). Such a tangent bifurcation 
in general produces a pair of stable and unstable solutions; here the solution with the smaller value of A, decreasing 
with K, should be unstable. Figure 7(c) shows that the unstable solution becomes null (A = 0) at VI w 1.257 as K 
approaches Kq « 3.515. Figure 7 also reveals the occurrence of the third transition at K c « 3.46, which is of the 
same nature as the second. 

The values of K c w 1.596 and K rj 3.515 can also be obtained analytically since they are given by the solutions 
of Eq. (|2|) in the limit A — > 0. In this limit, the right-hand side of the second equation vanishes, yielding fl — titt/t 
with n integer. The first one, which reduces to K~ x cos fir = (7r/2)g(f2), then gives 

K= -, y (23) 



where it has been noted that K > 0. Taking n — and n = 1 in Eq. (23), where g(u>) is given by the Gaussian 
distribution with unit variance, we obtain K = K c = w 1.596 (with Q = 0) and K = Kq = y^S/7re 2( - 7 ' 'l T ^> w 

3.515 (with Q = 2tt/t w 1.257) for r = 5, respectively. 

To examine how the stability changes at these bifurcations, we now consider a small perturbation from the incoherent 
state, for which the stationary distribution in Eq. (j^) is simply given by 1/2tt, and write 

P(lM) = -^+e>7foM), (24) 
where e 1. Upon substitution into Eq. (Q) and with A = Aie + 0(e 2 ), we obtain, to the lowest order in e, 

^ = -^ + ^A lC os(V>-0o), (25) 
at dip 2ir 

and seek solutions of the form 

» 7 (^,t) = c(t;w)e < * + c*(t;w)e- < * 1 (26) 
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where higher harmonics have been neglected. Equations ( |25| ) and fl26|), together with Eq. (|8j), lead to the amplitude 
equation for c(i;cD): 

dc ^^ = c (t, u) + ^-e lSlT r d£ug{£o + Sl)c(t-r, Q), (27) 

2 J-oo 

which in general possesses both discrete and continuous spectra. To find out the discrete spectrum, we put 

c(t]Q) = b{Cj)e xt , (28) 
where the eigenvalue A is independent of ui, and obtain the equation 

-{X-iQ)r K r duj iM - i ( 2 g) 

which has been examined for a Lorenzian distribution fl7|| . 

Here we investigate Eq. (|2^) for a Gaussian distribution. The stability of the incoherent state depends on whether 
all roots of Eq. (|29| ) possess negative real parts, i.e., Re A < 0. This is the case for K less than K s , where the 
incoherent state is neutrally stable (since the continuous spectrum is pure imaginary). Beyond K Sl there appears an 
eigenvalue with a positive real part, giving rise to instability. The value of K s can be computed from Eq. (E9) with 
Re A = imposed; this yields the coupled equations for K s and Im A 



cos 7 

v ° 

sta7T = 4 7e -^g_I_, (30) 

where 7 = £1 — Im A. Unlike the system without delay, Eq. ([5(j|) has an infinite number of solutions, among which the 
lowest value of K s should be taken. It is obvious that 7 = is the desired solution (regardless of time delay), leading 
to the lowest value K s — « 1.596. Note also that this value of K s coincides exactly with that of K c , implying 

that the incoherent state becomes unstable simultaneously with the appearance of the (stable) coherent state with 

n = 0. 

These results reveal that the order parameter exhibits a supercritical bifurcation at K = K c along the leftmost 
curve (£1 = 0) in Fig. 3. Namely, the emergence of a nontrivial solution (A > 0) is accompanied by the loss of 
stability of the null solution (A = 0) at K c (fl = 0) « 1.596. For the rest (£1 > 0), on the other hand, the unstable 
solution, generated together with the stable one by a tangent bifurcation at K c (tt), decreases as K is raised further 
and vanishes to zero at a larger value, K — Kq(U). For example, the unstable solution for O « 1.09, emerging at 
K w 1.97, decreases to zero at K 3.515 [see Fig. 7(c)]. It is thus concluded that for £1 > the bifurcation at Kq is 
subcritical: Between K c and Kq there exist an unstable coherent state in addition to the stable coherent states (and 
the incoherent one) although the unstable states have not been displayed in Fig. 3. 

The general features of the synchronization behavior obtained here are similar to those in Ref. Jl7| , and it is 
thus concluded that the difference in the distribution of natural frequencies does not change results qualitatively. 
On the other hand, we have examined additional interesting phenomena such as frequency suppression and details of 



multistability. In particular, unlike in the system with too ^ considered mostly in Ref. 17 1, here the phase boundaries 
with different values of the synchronization frequency do not intersect with each other on the K—t plane. [Compare 
Fig. 2(a) and (b).] Accordingly, the system with luq = does not undergo a discontinuous transition directly from the 
stable incoherent state and the coherent one with a nonzero synchronization frequency, and the associated hysteresis 
may not be observed. Further, in order to confirm these results, we have also performed numerical simulations, the 
results of which are presented in the next section. 



IV. NUMERICAL SIMULATIONS 

We have studied directly the equations of motion given by Eq. ([!]) via numerical simulations. The globally coupled 
system of size N = 5000, where natural frequencies are distributed according to the Gaussian distribution with unit 
variance, has been considered, and the Euler method with discrete time steps of St = 0.01 has been employed. At 
each run, we have discarded the first 10 5 time steps per oscillator to eliminate transient effects and taken the next 10 5 
time steps per oscillator to investigate synchronized solutions. Finally, independent runs with 30 different realizations 
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of the natural frequency distribution and initial conditions have been performed, over which the averages have been 
taken. In the simulations, the synchronization frequency is given by the average phase speed, i.e., the average rate 
of the phase change, and the obtained data at the coupling strength K = 10 are represented by crosses in Fig. 8(a). 
Note that both the incoherent state and the coherent one are found to be stable at the same value of r, indicating 
multistability; frequency suppression with increasing delay is also manifested. For comparison, the results shown 
in Fig. 1(a), obtained from Eq. (|l^), are also displayed, and perfect agreement is observed. Notice here that the 
basin of attraction shrinks rapidly with the synchronization frequency f2, which makes it quite difficult in numerical 
simulations to find the coherent-state solutions with large values of O. Figure 8(b) shows the behavior of the order 
parameter A as a function of the coupling strength K for r = (plus signs) and r = 5 (crosses). In both cases the 
system displays a continuous transition to the coherent state (with zero synchronization frequency). Slight suppression 
of synchronization by time delay can be observed. The error bars have been estimated by the standard deviation 
and the lines are guides to the eye. To make comparison of the analytical results obtained from Eq. (|l^) and the 
simulation results, we have also included in Fig. 8(b) the analytical results for r = 5, which are represented by the 
solid line. Good overall agreement between the two can be observed. 



We have studied analytically and numerically the collective synchronization phenomena in a set of globally coupled 
oscillators with time retarded interaction. In order to understand the effects of time delay on the synchronization, 
we have derived the self-consistency equations for the order parameter, which describe synchronization in the system. 
The detailed analysis of the self-consistency equations has revealed a multitude of coherent states with nonzero 
synchronization frequencies, each separated from the incoherent state by a discontinuous transition. At the critical 
coupling strength, the system exhibits the usual continuous transition from the incoherent state to the coherent 
one, displaying collective synchronization with zero synchronization frequency. As the coupling strength is increased 
further, coherent states with larger synchronization frequencies have also been shown to appear via discontinuous 
transitions from the incoherent state. Thus a multitude of coherent states with different synchronization frequencies 
have been found to coexist in the appropriate regions, leading to multistability. The synchronization frequency of the 
oscillators in a coherent state has been observed to decrease with the delay time. 

To confirm the analytical results, we have also performed numerical simulations, the results of which indeed display 
multistability and suppression of the synchronization frequency. For detailed comparison, however, one should search 
the solution space extensively, with varying initial conditions, to obtain solutions with various values of the synchro- 
nization frequency. This requires more extensive simulations, which is left for future study. Finally, one may also 
include stochastic noise in the system and study its effects on synchronization behavior. In particular, the interplay 
between the external driving and noise poses the possibility of stochastic resonance , and it is of interest to examine 
how the collective synchronization together with the time delay affects the possible resonance phenomena. 
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APPENDIX A 



In the case 51 <§; KA (<C 1), we approximate the integral appearing in Eq. (HO 
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and expand Eq. ( |l6| ) to the order (if A) 3 . This yields 
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where Eqs. ( p.9[ ) and the Gaussian distribution g(u>) have been used. After a tedious calculation, we obtain from 
Eq. (O) 
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APPENDIX B 

In the case f2 3> KA (^C 1), we approximate the integral in Eq. ( ]ilf ) as follows: 



9 



J dxg(fl+KAx) (x - \/ x 2 - lj + j dxg(Sl+KAx) (x + \J x 2 - lj 
= J dx [g(Q+KAx) - g(Q-KAx)} (x - \/ x 2 - I 
w - ^ dx g(n-KAx) (x - \/x 2 - lj , 

which, upon expansion to the order (AA) 3 , gives Eq. ( |l6| ) in the form 
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with the error function 4>(y) = (2/^/tt) J q dz e z . After a tedious calculation, we obtain from Eq. (B.2) 
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FIG. 1. Synchronization frequency Q (a) as functions of the delay time r at the coupling strength K — 10, where frequency 
suppression and multistability can be observed; (b) as functions of K at r = 5. For given K, the largest synchronization 
frequency belonging to the highest stair gives the most stable solution, but the corresponding basin of attraction is small. 
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FIG. 2. Phase diagram on the K-r plane, displaying boundaries between the incoherent and coherent states. Whenever 
each boundary in (a) is crossed from below, a new additional coherent state with a larger synchronization frequency emerges. 
Below the lowest boundary, which is the straight line K = K c m 1.596, only the incoherent state exists. In (b) the boundaries 
are redrawn, with the horizontal axis rescaled: f in (b) corresponds to Qt/(Q + 3) in (a). Here displayed is only the envelope 
consisting of the curve segments with lowest values of K at given f , which describes the phase boundary separating the 
incoherent state in the system with luq — 3. 



14 











4 

K 



8 



FIG. 3. Order parameter A as a function of the coupling strength K at the time delay r = 5. Each line describes the 
transition for each synchronization frequency, and the critical coupling strength K c is shown to increase for large fi. The 
leftmost curve, which corresponds to Q = 0, gives K c ~ 1.596; the next one, corresponding to Q ~ 1.09, gives K c « 1.97. The 
numerical results displayed by the leftmost curve agree well with the analytical ones. 
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FIG. 4. Graphical solutions of Eq. (§o|), displaying /(A) = {axK - 1)A + bi(is:A) 2 + ci(A"A) 3 versus A, for h = with (a) 
K — 1.594 (< K c ), (b) K — 1.596 (= K c ), and (c) K = 1.597 (> K c ). The negative solution appearing in (c) is unphysical. 
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FIG. 5. Synchronization frequency f2 versus the order parameter A obtained from Eq. © for r = 5 at (a) A' = 1.58, (b) 
K — 1.596, (c) K — 1.60, (d) K = 1.75. Solid and dotted lines represent solutions of the first equation and the second equation 
of Eq. (|l6|). A continuous transition for Q = can be observed at K — 1.596. 
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FIG. 6. Synchronization frequency SI versus the order parameter A obtained from Eq. (|l|) for t = 5 at (a) A' = 1.89, (b) 
K = 1.97, (c) K = 2.10, (d) K = 3.00. A discontinuous transition for Q « 1.09 can be observed at K = 1.97. 
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FIG. 7. Synchronization frequency f2 versus the order parameter A obtained from Eq. (|l|) for r = 5 at (a) if 
if = 3.46, (c) K = 3.515, (d) K — 4.00. A discontinuous transition for SI ~ 2.25 can be observed at K — 3.46. 
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FIG. 8. Results of numerical simulations on 5000 coupled oscillators: (a) Synchronization frequency versus the delay time. 
Crosses are results of numerical simulations and solid lines represent the solutions of Eq. (^6|), displaying perfect agreement with 
each other, (b) Order parameter versus the coupling strength, displaying continuous transitions (with zero synchronization 
frequency) for r = (plus signs) and r = 5 (crosses). The solid line represents analytic results for r = 5, displaying reasonable 
agreement with the numerical ones. Slight suppression of synchronization by time delay can be observed near the transition 
region. The size of the error bars estimated by the standard deviation is about that of the symbols and lines are merely guides 
to the eye. 
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